/*
 * Copyright 2010-2012 Susanta Tewari. <freecode4susant@users.sourceforge.net>
 *
 * This program is free software: you can redistribute it and/or modify
 * it under the terms of the GNU General Public License as published by
 * the Free Software Foundation, either version 3 of the License, or
 * (at your option) any later version.
 *
 * This program is distributed in the hope that it will be useful,
 * but WITHOUT ANY WARRANTY; without even the implied warranty of
 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 * GNU General Public License for more details.
 *
 * You should have received a copy of the GNU General Public License
 * along with this program.  If not, see <http://www.gnu.org/licenses/>.
 */

package bd.org.apache.commons.math.linear;

import bd.org.apache.commons.math.exception.DimensionMismatchException;
import bd.org.apache.commons.math.exception.NonSquareMatrixException;
import bd.org.apache.commons.math.util.FastMath;

/**
 * Calculates the LUP-decomposition of a square matrix.
 * <p>The LUP-decomposition of a matrix A consists of three matrices L, U and
 * P that satisfy: P&times;A = L&times;U. L is lower triangular (with unit
 * diagonal terms), U is upper triangular and P is a permutation matrix. All
 * matrices are m&times;m.</p>
 * <p>As shown by the presence of the P matrix, this decomposition is
 * implemented using partial pivoting.</p>
 * <p>This class is based on the class with similar name from the
 * <a href="http://math.nist.gov/javanumerics/jama/">JAMA</a> library.</p>
 * <ul>
 * <li>a {@link #getP() getP} method has been added,</li>
 * <li>the {@code det} method has been renamed as {@link #getDeterminant()
 * getDeterminant},</li>
 * <li>the {@code getDoublePivot} method has been removed (but the int based
 * {@link #getPivot() getPivot} method has been kept),</li>
 * <li>the {@code solve} and {@code isNonSingular} methods have been replaced
 * by a {@link #getSolver() getSolver} method and the equivalent methods
 * provided by the returned {@link DecompositionSolver}.</li>
 * </ul>
 *
 * @version $Id: LUDecomposition.java 1244107 2012-02-14 16:17:55Z erans $
 * @see <a href="http://mathworld.wolfram.com/LUDecomposition.html">MathWorld</a>
 * @see <a href="http://en.wikipedia.org/wiki/LU_decomposition">Wikipedia</a>
 * @since 2.0 (changed to concrete class in 3.0)
 */
public class LUDecomposition {

    /**
     * Default bound to determine effective singularity in LU decomposition.
     */
    private static final double DEFAULT_TOO_SMALL = 1e-11;

    /**
     * Entries of LU decomposition.
     */
    private final double[][] lu;

    /**
     * Pivot permutation associated with LU decomposition.
     */
    private final int[] pivot;

    /**
     * Parity of the permutation associated with the LU decomposition.
     */
    private boolean even;

    /**
     * Singularity indicator.
     */
    private boolean singular;

    /**
     * Cached value of L.
     */
    private RealMatrix cachedL;

    /**
     * Cached value of U.
     */
    private RealMatrix cachedU;

    /**
     * Cached value of P.
     */
    private RealMatrix cachedP;

    /**
     * Calculates the LU-decomposition of the given matrix.
     * This constructor uses 1e-11 as default value for the singularity
     * threshold.
     *
     * @param matrix Matrix to decompose.
     */
    public LUDecomposition(RealMatrix matrix) {

        this(matrix, DEFAULT_TOO_SMALL);
    }

    /**
     * Calculates the LU-decomposition of the given matrix.
     *
     * @param matrix The matrix to decompose.
     * @param singularityThreshold threshold (based on partial row norm)
     * under which a matrix is considered singular
     */
    public LUDecomposition(RealMatrix matrix, double singularityThreshold) {

        if (!matrix.isSquare()) {

            throw new NonSquareMatrixException(matrix.getRowDimension(),
                                               matrix.getColumnDimension());
        }

        final int m = matrix.getColumnDimension();

        lu      = matrix.getData();
        pivot   = new int[m];
        cachedL = null;
        cachedU = null;
        cachedP = null;

        // Initialize permutation array and parity
        for (int row = 0; row < m; row++) {

            pivot[row] = row;
        }

        even     = true;
        singular = false;

        // Loop over columns
        for (int col = 0; col < m; col++) {

            // upper
            for (int row = 0; row < col; row++) {

                final double[] luRow = lu[row];
                double         sum   = luRow[col];

                for (int i = 0; i < row; i++) {

                    sum -= luRow[i] * lu[i][col];
                }

                luRow[col] = sum;
            }

            // lower
            int    max     = col;    // permutation row
            double largest = Double.NEGATIVE_INFINITY;

            for (int row = col; row < m; row++) {

                final double[] luRow = lu[row];
                double         sum   = luRow[col];

                for (int i = 0; i < col; i++) {

                    sum -= luRow[i] * lu[i][col];
                }

                luRow[col] = sum;

                // maintain best permutation choice
                if (FastMath.abs(sum) > largest) {

                    largest = FastMath.abs(sum);
                    max     = row;
                }
            }

            // Singularity check
            if (FastMath.abs(lu[max][col]) < singularityThreshold) {

                singular = true;

                return;
            }

            // Pivot if necessary
            if (max != col) {

                double         tmp   = 0;
                final double[] luMax = lu[max];
                final double[] luCol = lu[col];

                for (int i = 0; i < m; i++) {

                    tmp      = luMax[i];
                    luMax[i] = luCol[i];
                    luCol[i] = tmp;
                }

                int temp = pivot[max];

                pivot[max] = pivot[col];
                pivot[col] = temp;
                even       = !even;
            }

            // Divide the lower elements by the "winning" diagonal elt.
            final double luDiag = lu[col][col];

            for (int row = col + 1; row < m; row++) {

                lu[row][col] /= luDiag;
            }
        }
    }

    /**
     * Returns the matrix L of the decomposition.
     * <p>L is a lower-triangular matrix</p>
     *
     * @return the L matrix (or null if decomposed matrix is singular)
     */
    public RealMatrix getL() {

        if ((cachedL == null) &&!singular) {

            final int m = pivot.length;

            cachedL = MatrixUtils.createRealMatrix(m, m);

            for (int i = 0; i < m; ++i) {

                final double[] luI = lu[i];

                for (int j = 0; j < i; ++j) {

                    cachedL.setEntry(i, j, luI[j]);
                }

                cachedL.setEntry(i, i, 1.0);
            }
        }

        return cachedL;
    }


    /**
     * Returns the matrix U of the decomposition.
     * <p>U is an upper-triangular matrix</p>
     *
     * @return the U matrix (or null if decomposed matrix is singular)
     */
    public RealMatrix getU() {

        if ((cachedU == null) &&!singular) {

            final int m = pivot.length;

            cachedU = MatrixUtils.createRealMatrix(m, m);

            for (int i = 0; i < m; ++i) {

                final double[] luI = lu[i];

                for (int j = i; j < m; ++j) {

                    cachedU.setEntry(i, j, luI[j]);
                }
            }
        }

        return cachedU;
    }


    /**
     * Returns the P rows permutation matrix.
     * <p>P is a sparse matrix with exactly one element set to 1.0 in
     * each row and each column, all other elements being set to 0.0.</p>
     * <p>The positions of the 1 elements are given by the {@link #getPivot()
     * pivot permutation vector}.</p>
     *
     * @return the P rows permutation matrix (or null if decomposed matrix is singular)
     * @see #getPivot()
     */
    public RealMatrix getP() {

        if ((cachedP == null) &&!singular) {

            final int m = pivot.length;

            cachedP = MatrixUtils.createRealMatrix(m, m);

            for (int i = 0; i < m; ++i) {

                cachedP.setEntry(i, pivot[i], 1.0);
            }
        }

        return cachedP;
    }


    /**
     * Returns the pivot permutation vector.
     *
     * @return the pivot permutation vector
     * @see #getP()
     */
    public int[] getPivot() {

        return pivot.clone();
    }


    /**
     * Return the determinant of the matrix
     *
     * @return determinant of the matrix
     */
    public double getDeterminant() {

        if (singular) {
            return 0;
        } else {

            final int m           = pivot.length;
            double    determinant = even
                                    ? 1
                                    : -1;

            for (int i = 0; i < m; i++) {

                determinant *= lu[i][i];
            }

            return determinant;
        }
    }


    /**
     * Get a solver for finding the A &times; X = B solution in exact linear
     * sense.
     *
     * @return a solver
     */
    public DecompositionSolver getSolver() {

        return new Solver(lu, pivot, singular);
    }


    /**
     * Specialized solver.
     */
    private static class Solver implements DecompositionSolver {

        /**
         * Entries of LU decomposition.
         */
        private final double[][] lu;

        /**
         * Pivot permutation associated with LU decomposition.
         */
        private final int[] pivot;

        /**
         * Singularity indicator.
         */
        private final boolean singular;

        /**
         * Build a solver from decomposed matrix.
         *
         * @param lu entries of LU decomposition
         * @param pivot pivot permutation associated with LU decomposition
         * @param singular singularity indicator
         */
        private Solver(final double[][] lu, final int[] pivot, final boolean singular) {

            this.lu       = lu;
            this.pivot    = pivot;
            this.singular = singular;
        }

        /**
         * {@inheritDoc}
         */
        @Override
        public boolean isNonSingular() {

            return !singular;
        }


        /**
         * {@inheritDoc}
         */
        @Override
        public RealVector solve(RealVector b) {

            final int m = pivot.length;

            if (b.getDimension() != m) {
                throw new DimensionMismatchException(b.getDimension(), m);
            }

            if (singular) {
                throw new SingularMatrixException();
            }

            final double[] bp = new double[m];

            // Apply permutations to b
            for (int row = 0; row < m; row++) {

                bp[row] = b.getEntry(pivot[row]);
            }

            // Solve LY = b
            for (int col = 0; col < m; col++) {

                final double bpCol = bp[col];

                for (int i = col + 1; i < m; i++) {

                    bp[i] -= bpCol * lu[i][col];
                }
            }

            // Solve UX = Y
            for (int col = m - 1; col >= 0; col--) {

                bp[col] /= lu[col][col];

                final double bpCol = bp[col];

                for (int i = 0; i < col; i++) {

                    bp[i] -= bpCol * lu[i][col];
                }
            }

            return new ArrayRealVector(bp, false);
        }


        /**
         * {@inheritDoc}
         */
        @Override
        public RealMatrix solve(RealMatrix b) {

            final int m = pivot.length;

            if (b.getRowDimension() != m) {
                throw new DimensionMismatchException(b.getRowDimension(), m);
            }

            if (singular) {
                throw new SingularMatrixException();
            }

            final int nColB = b.getColumnDimension();

            // Apply permutations to b
            final double[][] bp = new double[m][nColB];

            for (int row = 0; row < m; row++) {

                final double[] bpRow = bp[row];
                final int      pRow  = pivot[row];

                for (int col = 0; col < nColB; col++) {

                    bpRow[col] = b.getEntry(pRow, col);
                }
            }

            // Solve LY = b
            for (int col = 0; col < m; col++) {

                final double[] bpCol = bp[col];

                for (int i = col + 1; i < m; i++) {

                    final double[] bpI    = bp[i];
                    final double   luICol = lu[i][col];

                    for (int j = 0; j < nColB; j++) {

                        bpI[j] -= bpCol[j] * luICol;
                    }
                }
            }

            // Solve UX = Y
            for (int col = m - 1; col >= 0; col--) {

                final double[] bpCol  = bp[col];
                final double   luDiag = lu[col][col];

                for (int j = 0; j < nColB; j++) {

                    bpCol[j] /= luDiag;
                }

                for (int i = 0; i < col; i++) {

                    final double[] bpI    = bp[i];
                    final double   luICol = lu[i][col];

                    for (int j = 0; j < nColB; j++) {

                        bpI[j] -= bpCol[j] * luICol;
                    }
                }
            }

            return new Array2DRowRealMatrix(bp, false);
        }


        /**
         * {@inheritDoc}
         */
        @Override
        public RealMatrix getInverse() {

            return solve(MatrixUtils.createRealIdentityMatrix(pivot.length));
        }
    }
}
